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There is good evidence from N-body simulations that the velocity distribution in the outer parts of halos is 
radially anisotropic, with the kinetic energy in the radial direction roughly equal to the sum of that in the two 
tangential directions. We provide a simple algorithm to generate such cosmologically important distribution 
functions. Introducing r^iE), the radius of the largest orbit of a particle with energy E, we show how to write 
down almost trivially a distribution function of the form /(£, L) = L^^gir^) for any spherical model - including 
the 'universal' halo density law (Navarro-Frenk- White profile). We in addition give the generic form of the 
distribution function for any model with a local density power-law index a and anisotropy parameter p and 
provide limiting forms appropriate for the central parts and envelopes of dark matter halos. From those, we 
argue that, regardless of the anisotropy, the density falloff at large radii must evolve to p ~ r^* or steeper 
ultimately. 

PACS numbers: 95.35. +d, 98.62.Gq 



I. INTRODUCTION 

N-body experiments now can reliably follow the collapse 
and violent relaxation of dark matter halos from initial con- 
ditions. This has led to the discovery of regularities in the 
phase space distribution of dark matter [e.g., 1], even though 
the final state is not completely independent of initial con- 
ditions. This is important because it suggests that there is a 
generic functional form for the distribution function (DF) that 
describes the physics of violent relaxation, albeit with some 
cosmic scatter |2]. 

For example, Hansen and Moore ||3|, see also01 have found 
that the density power index is correlated with the anisotropy 
parameter j6 = 1 - (v^)/(2(v?)) f3\. Here, <v^) and <v^) are the 
radial and the tangential velocity second moments. For a wide 
range of cosmological simulations, they demonstrate that the 
dark matter follows the equation of state/? ^ I- 1.15(1 -a/6) 
where a is the density power index (i.e., p ~ r In the very 
center, dark matter halos are roughly isotropic (J3 ~ Q) with 
ax 1 . In the outer parts, violent relaxation produces a density 
profile that asymptotically becomes p ~ r'^ |l6] orp ~ r"^ 01, 
for which the anisotropy parameter/? ^ 0.5 accordingly. 

If violent relaxation proceeded to completion, then equipar- 
tition would enforce equal kinetic energy in each direction and 
the velocity distribution would be isotropic |8]. This appears 
to be the case only at the centers of numerical simulations. 
Particles with large apocenters respond only weakly to the 
fluctuating gravitational field. Throughout most of the halo, 
this gives rise to an end point for which the kinetic energy 
in the radial direction is roughly equal to the sum of that in 
the two tangential directions. This seems to be supported not 
only by the numerical simulations but also by the observation 
of stars in elliptical galaxies ||9|], whose kinematics is also gov- 
erned by the collisionless Boltzmann equation with the gravi- 
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tational potential. The purpose of this paper is to give the DF 
of the dark matter which has this property. 

There has been much work on isotropic DFs [see 5] of grav- 
itating systems. These are fine for the inner parts. On the other 
hand, there has been much less work on DFs suitable for the 
radially anisotropic outer parts of the dark matter halos. In 
particular, a number of the suggestions in the literature for 
anisotropic DFs [e.g., are unsuitable, as they 

yield overwhelming radial anisotropy (J3 ^ 1) in the outer 
parts, which is inconsistent with the simulations. While there 
exist some suggestions on the form of anisotropic DFs with a 
more flexible behavior of /? [e.g., , recovering such DFs 

for most density profiles is often analytically intractable 
but see 12 for a special case]. 



II. DISTRIBUTION FUNCTIONS WITH /3 = 1 /2 

The widely used ansatz for a DF of a spherical system with 
constant anisotropy (parameterized by /3) is 



f{E,L)^L-^PfE{E) 



(1) 



where E - t// - /2 is the binding energy per unit mass, L = 
rvj is the specific angular momentum, and t// is the relative 
potential. Integration of the DF over the velocity gives 
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The unknown function /^(E) then can be recovered from the 
integral inversion formula 11131 Usll : 



fE(E) = 
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where h = r^p is expressed as a function of ifr, and n - [(3 /2- 
/3)i and A = 3/2- p-n are the integer floor and the fractional 
part of 3/2 - /?. This includes Eddington's formula ifl^ for 
the isotropic DF as a special case (J3 = 0). The expression for 
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the differential energy distribution (DED) reduces to Fcf.. 11311 



dM 
'dE 



fE(E) 
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Here, is the radius of the largest orbit of a particle with 
energy E, that is to say, iffirg) = E. 

If yS is a half-integer constant (i.e., /3 - 1/2, -1/2, and so 
on), the expression for DF further reduces to 



f(E,L). 
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This involves only differentiations, as first noted by Cudde- 
ford lilQl- For the simplest case of yS = 1/2, by utilizing the 
parameter r^, the expressions for the DF and the DED can be 
simply written down as 



f(E,L) = 
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where Mr is the enclosed mass within the sphere of radius of 
r. So, the (3—1/2 case, which is desirable from the point of 
view of the N-body simulations, is also very attractive mathe- 
matically. That is, both the DF and DED can be found almost 
trivially from the potential-density pair 



III. COSMOLOGICAL HALO MODELS 



A. Generalized NFW Profiles 



Let us consider a family of centrally cusped density profiles 
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AttG r(r + af-'^ 
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where the parameter > 2 is the asymptotic density power 
index at large radii, and x/zq is the depth of the central potential 
well. Near the center, the density for every member of this 
family is always cusped as r This reduces to the Hernquist 
I2OI model for b - A- whereas it becomes the 'universal' halo 
density law or so-called Navarro-Frenk- White (NFW) profile 
iP] if = 3. The system has an infinite mass for 2 < b < 3. 
On the other hand, if b > 3, the finite total mass is given by 
GMoo - if/Qa/ib - 3). The corresponding potential is 



^^Ln,_,(^), 
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where, to reduce notational clutter, we have used the 
logarithm' function [2 1.1 defined to be 



Ln„(x) 



^ J^'^_|(x'-«-l)/(l-^) 
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Here, we note also a property of q'-logarithm function, namely 
Ln^(jic"') = -Ln2-,(x). Its inverse is the '^-exponential' func- 
tion 

Ex^(.x:) ^\ja. '(x)^\ 

' I exp(jc) q - 1 

The DF of the form of Eq. Q can be found using (with 
G - tpi) - a - I) 



r^Pp-. 
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and Eqs. (0 or However, for y6 = 1/2, the formulas (|6j 
enable us to write down the DF for all the family using 
[here, r^Zi = Lni_2(l + >'e)\ see AppendixIXIfor b = 7/3, 5/2, 
8/3,7/2, 4, or 5]. We find that 
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Here, E ranges in the interval [0, 1] because < E < tfr < 1. 
Note that — > as £ — > 1 and — > 00 as £ — » 0. For all 
members of the family with b>2,we find that the DF is non- 
negative for all accessible phase space volume. The behavior 
of the energy part of the DF (Eq.|8} and the DED (Eq.|9|l of 
this family are shown in Fig.^for several values of b. 

Lokas and Mamon [22} derived analytical expressions for 
various physical properties of the NFW profile including the 
profiles of the velocity dispersions and the kinetic and po- 
tential energy assuming isotropy, constant anisotropy, or an 
Osipkovl 10, ll]-MeiTitt| 12| type DF. However, they fell short 
of deriving any explicit DF with the exception of the isotropic 
DF, for which they gave the result of the numerical integration 
of the Eddington's formula. In fact, (the numerical integra- 
tions of) the isotropic as well as the Osipkov-Merritt DFs and 
DEDs for the NFW profile have been investigated in detail 
by Widrow |23]. We nevertheless again note that the isotropic 
DF is suitable only for the very inner parts of dark matter halos 
and that the Osipkov-Merritt type DF is unsuitable in general, 
as it yields overwhelming radial anisotropy (/? — » 1) in the 
outer parts. 

It is relatively straightforward to find the asymptotic be- 
havior of the DF and the DED near £ = 1 by means of 
Taylor series expansion of at £ = 1. We find that 
f{E = 1,L) = 2{b- 1)/[(27t)^L] tends to a constant while 
dM/dE ~ (1 - £')2 -> as £ -> 1. On the other hand, the 
asymptotic behavior near E - for b ^ 3 can be derived from 
(3 - b)rEE X rlr'' - 1 for r£ » 1 . Then, we have 



rE 



^-\l(h-2) 2<b<3 
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for < £ <K 1 . Consequently, the asymptotic forms of the DF 
and the DED are given by 
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FIG. 1 : The energy part of the distribution function (upper panel) and the differential energy distribution (lower panel) of the Navan'o-Frenk-White family with 
the constant anisotropy parameter of /? = 1/2: dotted lines {b = 5/2), thin solid lines (h = 3; NFW profile), short-dashed lines (b = 7/2), think solid lines (b = 4; 
Hemquist model), long-dashed lines (b = 9/2), dot-dashed lines {b = 5). Note that the models are normalized to the common value of the depth of the central 
potential well, and, for b < 3, the total mass is divergent. 



SO that f(E - 0,L) - whereas the DED diverges as £ —> if 
b <4 and is finite otherwise. In particular, limE^o(dM/dE) - 
if /7 > 4, and dM/dE\E^Q = 3 if = 4. 

For the NFW profile (b = 3), the proper asymptotic form for 
the inversion of £ = r^' ln(l -i- r^) cannot be expressed using 
only elementary functions. Nevertheless, the continuous na- 
ture of the asymptotic behavior suggests that dM/dE ~ 
In fact, since E ~ r^' Inr^ and dM/dE ~ rElirzE - 1) ~ 
rE/\nrE as — > oo, this is indeed the right behavior for 
the NFW profile. We also note that it is possible to approxi- 
mate as rE ~ E-^Q.{E) and /e ~ E[D.(E)r^ where □(£) = 
-W-ii-E) = ln(£-i ln(£-i ln(£-' . . .))) ^ ln(£-i Infi-'). 
Here, W-i(x) is the real-valued second branch of Lambert W- 
function 12411 such that W-i(x) < -1 for -e ' < x < 0. By 
comparison, we have / ~ E^^^[Q.{E)]^^ and dM/dE ~ 
for the same {E — > 0)-asymptotic behaviors of the isotropic 
DF and DED of the NFW profile lEl . 

As a brief illustration of the application of our DF (Eq. 
we consider the direct detection of dark matter. Such exper- 
iments work by measuring the recoil energy of a nucleus in 
a low background laboratory detector that has undergone a 
collision with a dark matter particle. Although the deposited 
energy is tiny and the interactions are very rare, there are now 
many groups searching for this eff'ect worldwide [e.g.,|25l and 
references therein]. The detection rate depends on the masses 
TO^ and Ton of the dark matter particle and the target nucleus 
and the elastic scattering cross-section ctq between them. It 
also depends on the local dark matter density po and the speed 
distribution of the dark matter particles (in the rest frame of 
the target). 



The differential rate for detection per unit detector mass is 
given by [c.f.,|2^ 



dR 
dE 



e=e 2to^jU 



Cd^v 
J 1^^ 
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where 0(x) is the Heaviside unit step function, £ is the recoil 
energy,//"' = to,^' + m^^ is the reduced mass, and 

In addition, F(S) is the nuclear form factor, which is com- 
monly modeled at least for scalar interaction by L27lE3l 



F(fi) = exp - 



where Go is the nuclear coherence energy. Here, the DF is 
normalized to the local density of the dark matter; 



PO 



d\f(E,L). 



Note the integral in Eq. ilQ\ is over the velocity with respect to 
the detector on Earth. Of course, the Earth revolves around the 
Sun while the Sun moves with respect to the Galactic inertial 
frame (in which the net angular momentum of the dark matter 
halo vanishes). This produces an annual modulation in the 
signal, which the experiments hope to detect. 

The total event rate can be found by integrating over all 
detectable energies. For the sake of definiteness, we consider 
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FIG. 2: Upper panel: The differential rate (in units of events kg"' d"' keV"') 
for the case of 40 GeV dark matter particles impinging on a cryogenic detec- 
tor made of Ge. Results are shown for a standard isothermal sphere with 
Maxwellian DF (full line) and the NFW model with DF given by Eq. jS) 
(dashed line). Lower panel: The annual modulation signal (in units of events 
kg" ' d" ' ) for the isothermal sphere and NFW models. The upper curves show 
the variation in the total rate, the lower curves the variation in the low energy 
events (< 10 keV). 



a dark matter particle with only scalar interactions and with a 
mass m-^c- - 40 GeV and cross-section o-q = 4 x lO"-'^ cm"^. 
The detector is made of ^^Ge. The local halo density is taken 
aspoc^ =0.3 GeV cm^^ gl. 

We consider two models for the dark matter halo. The first 
is a standard isothermal sphere, with a flat rotation curve of 
amplitude vo - 220 km s The DF is a Maxwell-Boltzmann 
distribution |5|. This is a useful benchmark, as the model is 
widely used in dark matter studies [e.g., 26]. The second is a 
NFW model (Eq.0with b - 3) with a = 10 kpc, normalized 
to provide the assumed local halo density. The DF is given by 
Eq. (0. 

The results for the differential rate and the annual mod- 
ulation signal are shown in Fig. |2] The total rate is lower 
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FIG. 3: The energy part of the distribution function (upper panel) and the 
differential energy distribution (lower panel) of the y models with the constant 
anisotropy parameter ofyS = 1/2: solid lines (7 = 1; Hemquist model), dotted 
lines (y = 3/2), short-dashed lines (7 = 2; Jaffe model), long-dashed lines 
(y = 5/2). The models are normalized to the same total mass. For y > 2, 
the central potential well depth is infinite so that E lies in the range [0, 00), 
whereas < E < i/zq = (2 - yr\GM/a) fory < 2. 



by ~10% for the NFW model compared against the isother- 
mal sphere. But, the peak-to-peak amplitude of the modu- 
lation signal has increased from ~ 0.018 events kg"' d"' to 
~ 0.030 events kg"' d"'. This renders the dark matter par- 
ticle more detectable. However, this good news comes with 
a caveat. If the experiment is only sensitive to low energy 
events (£ < 10 keV) then the peak-to-peak variation in the 
modulation signal is actually smaller for the NFW model ver- 
sus the standard isothermal sphere. The main differences be- 
tween the two models is that the escape speed of dark mat- 
ter particles is finite for the NFW model, but infinite for the 
isothermal sphere. Therefore, the former produces more low 
energy events, while the latter provides a larger total number 
of events. 



B. The Gamma Spheres 

As a second example of our formulas (|6j, we consider the 
potential-density pair of the y model |30, 31] 



(3 - y)M 



GM^ lr + a\ 



An ry(r + af-y 



Here, the parameter y < 3 is the three-dimensional central 
density slope (i.e., the central density is cusped as r^^ if 7 > 
0), and M is the total mass. At large radii, the density falls off 
as r"'* for every member of the family. The central potential 
well depth is infinite if 2 < y < 3, whereas, if y < 2, the 
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potential is bounded as < < (2 - ■y)'^(GM/a). The famil 
contains the Hernquist |20| model (y - I) and the JafFe 
model (y = 2) as special cases. 

If we define y = r/(r + a), it is easy to write down (G 
M = a= 1) 



h = r^^p = 



s-rCi-y)' 



4;r 



Here, h can be written explicitly as a function of i/r using y — 
Exy_i (-)/'). Then, for [5 - 1/2, from Eq. (|5jl, we find the DF 
and the DED, 

f{E,L) = — ^[(4_^)3;£ + (^_l)] ; (11) 



(2;r)3L 
t/M _ 3-r 



-777 = -r^[(4-r)3'£ + (r-l)]3'£ 



(12) 



where yE - rzlil + rs) - Exy_i(-£'), which can be always 
expressible using elementary functions of E. Since < < 
1, the DF is everywhere non-negative only if 1 < y < 3 |33, 
c.f.,]. The behavior of the energy part of the DF (Eg. II 11 1 and 
the DED (Eq. I12> of the y models are shown in Fig. |3l for 
several values of y. 

As £■ — » (yE — » 1), we find that /^(-E) ~ E^ — > whereas 
the DED is always finite with the limiting value of 3(3 - y)/2. 
On the other hand, the asymptotic behavior as £ — > i/tq (yE — > 
0), where i/^q = oo for y > 2 or t//o - (2 - y)"' for y < 2, are 
found to be /e ~ y^' and (dM/dE) ~ yE where 



£■-1/(7-2) 



2 < y < 3 

r = 2 

[ (lAo - £)i/(2-r) 1 < y < 2 



except for the limiting case of y = 1 (the Hernquist model), 
for which /(£, L) = 3(471^)-^ L-^l - yE? = ^{'Xn^^E^L-^ 
and (dM/dE) = 3yE - 3(1 - E)^. Note that the asymptotic 
behavior of the DED at both limits are the same as those for 
the isotropic DF |30] for 1 < y < 3 despite the fact that the 
behavior of the DFs are rather distinct. 



IV. THE UNIVERSAL ASYMPTOTIC BEHAVIOR 

The asymptotic behaviors derived for the DF and DED in 
the preceding sections strongly suggest that they are simply 
determined by the anisotropy parameter and the density power 
index at the center and at the infinity. Even more interestingly, 
the DED appears to be completely determined (up to scale) by 
the density power index alone. Assuming that this is indeed 
the case, we predict asymptotic behaviors of generic DFs and 
DEDs by generalizing the method of Hjorth and Madsen 
while allowing for anisotropy by means of the ansatz Q. 

First, the self-consistent potential due to the asymptotic 
density profile of p ~ r " is 



r-(»-2) 



a > 3 
2 < ff < 3 



Here, since the case that o- > 3 is only allowed for the asymp- 
totic falloff at infinity, the potential should tend to the Kep- 
lerian finite mass limit. Next, if we assume /e to be roughly 
scale-free with /e ~ E" (for a > 2) or /£ ~ (i/^o - E)" (for 
a < 2), we find from Eq. (|3 

|^-2/?^«-/J+3/2 



In order for this to be p 
tial, we should have 



r " with the self-consistent poten- 



L'-^f{E,L) 



£/s-(i/2)+p 2 < c < 3 

exp[2(l -j8)£] a = 2 

(lAo - £)/'-(i/2)-(-/') 2/3 < a < 2 



(13) 



where p - 2(1 -/3)/{a-2). Here, a and /3 are the limiting val- 
ues at the center (for behavior near E ~ t//o) or the asymptotic 
value toward infinity (for behavior near E ~ 0). As for the 
DED, by changing the integration variable to ifr/E (for a > 2) 
or (iffo -iff)/ (iffo - E) (for a < 2) in Eq. ©, and combined with 
Eqs. O, 

we find that. 



dM 
'dE 



(E"-^ a>3 

^-\l(a-2) 2<a<3 

exp(-£') a - 2 

[{ipQ-Efl'-^-"'^ 2/3<a<2 



(14) 



(-Ao - -A) 



.,2-ff 
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It is independent of the anisotropy. We note that the coefficient 
for the leading order term of changes its sign aia -2/} so 
that the result is invalid at the limit a -2/5 and that the DF is 
unphysical for a < 2/3 ll33ll . 

Here, since or < 3 at the center, \\mE^^„ dM/dE - (where 
i/ro = ooif2<af<3) for all physical values of a. On the other 
hand, the behavior of dM/dE near E - Q for a finite mass sys- 
tem (i.e., a > 3>) implies that dM/dE diverges for a <A while 
it is finite for a > 4 (in particular, dM/dE — > if o- > 4). 
It has been argued before that violent relaxation produces an 
r~'^ density falloff at large radii jfl . This may be inferred from 
the generic behavior of dM/dE at the asymptotic limit £ — » 
(note that the behavior of dM/dE near £ = is dominated 
by the particles at large radii). That is, the loss of loosely 
bound particles at large radii due to velocity perturbations is 
much more significant if the initial density falloff is shallower 
than r"^ while it becomes rather insignificant once the den- 
sity falloff" gets steeper than r"^. Therefore, any perturbation 
drives systems with initially shallower density falloff to settle 
toward the r falloff or a slightly steeper slope. Note that this 
argument is completely independent of the anisotropy since 
the asymptotic form of dM/dE is also independent of (3. 

If one considers the case that the gravitational field is dom- 
inated not by the self-consistent potential but by the Keplerian 
field due to the central black hole, it is straightforward to see 
that the same argument leads to the asymptotic behavior of 
fE ~ E"-l^-^l^ and dM/dE ~ E"-'* for all allowed values of 
Of > /3 + 1 /2 except also for the limiting case a - [i + 1/2, for 
which the otherwise leading term identically vanishes. 
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V. CONCLUSIONS 

Except close to the very center, dark matter halos have dis- 
tribution functions (DFs) that are radially anisotropic with an 
anisotropy parameter p ^ 0.5. Constant anisotropy distribu- 
tion functions with yS = 1 /2 are very simple to construct - far 
simpler than Eddington's awkward Abel transformation pair 
for the isotropic model. This paper provides simple inversion 
formulas (|6|l, which enables such DFs and the corresponding 
differential energy distributions to be built for any spherical 
model, provided that its potential and density are known. 

We have also given the generic form of the DF for any 
spherical model with a local density power index a (p ~ r ") 
and anisotropy parameter y6. In the central parts of simulated 
halos where a ^ \ and a; 0, we find f(E) ~ (i^o - EY^^^ and 
dM/dE ~ (i^o - E). In the envelopes of simulated halos with 
y6 ~ 1/2, we have DFs, which range from 

IE dM I 

f(E, L) ; 

L A(E) dE E 

if a 3 [where A(£') is some function that diverges or van- 
ishes no faster than finite power of the logarithm as E — > 0], 



to 

£■2 dM n 
f(E, L) ; E° 

\f a K A. Finally, we have argued that falloff of the density at 
large radii must have evolved to p ~ r"** or steeper in the long 
run. This argument is independent of the velocity anisotropy. 

Using our DF for the NEW profile, we calculated the direct 
detection rate for a cosmological halo model with a radially 
anisotropic DF. We showed that the annual modulation signal 
is larger in radially anisotropic (j6 = 1 /2) cosmological halo 
models than in isotropic isothermal spheres. This may be wel- 
come good news for dark matter experimentalists. 



APPENDIX A 

Ifb^ 113, 5/2, 8/3, 7/2, 4, or 5, then the integral of motion 
rE for the generalized NEW profiles is expressible analytically 
in terms of the energy by solving a quadratic or linear equa- 
tion, viz {G - 4i(i - a - \) 



4(1 -£) 

9(3 - 4£2) + 33/2(3 _ 2£)3/2 Vl +2E 
A-E-E^l'^ iS + E 



rE = 



2E 



for b^- : 
2 



for b - — : 
3 



for b 



l-E 



rE = 



for ^7 = 4 ; 



rE 



rE 



33/2 yf4rE 9 

2£3/2 

1-4-E+ Vl +8E 
4E 



— for b - — : 

2E 3 



for = 5. 



For all other values of b, the function rE{E) is straightfor- 
ward to construct numerically. 
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